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In  the  theory  of  blast  vaves,  extensive  numerical  calculations  have 
become  a  standard  tool.  As  exen^jlified  by  the  collaborative  work  of  Boyer, 
Brode,  Glass  and  Hall  (1958) >  such  solutions  cam  provide  valuable  inter¬ 
pretation  for  e:i5>eri!aental  programs.  Conversely,  the  application  of  such 
theoretical  solutions  to  laboratory  or  field  experience  provides  a  desir¬ 
able  check  on  the  limitations  and  physical  approximati^ 3  of  the  theory. 
Kiese  factors  have  pron5)ted  the  present  paper,  in  which  £he  general  steps 


in  proceeding  from  physical  assumptions  to  numerical  answers  are  outlined 
for  solutions  by  a  Lagrangean  coordinate,  artificial  viscosity  method. 
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I.  imODUCTIOH 

During  the  past  tw  deceuies  there  has  been  an  increasing  ability  as 
well  as  a  continuing  interest  in  describing  the  hydrodynamic  consequences 
of  explosions.  The  history  of  shock  hydrodynamics  from  explosive  sources 
and  the  motivations  for  investigations  in  this  field  mke  an  interesting 
subject  for  discussion,  but  for  the  purposes  of  this  jwiper  attention  vill 
be  restricted  to  those  developments  vhich  have  led  to  the  current  ability 
to  carry  out  numerical  computations  describing  in  detail  the  hydrodynamic 
par’ameters  associated  with  explosions. 

After  the  fundamental  physics  governing  the  'Uieniiodynamlcs  and  hydro¬ 
dynamics  of  explosions  became  evident,  a  nvimber  of  investigators  were 
successiUl  in  formulating  analytical  approximations  for  e3q)losion  phenomena. 
One  type  of  solution,  which  has  been  variously  expressed  and  derived,  is 
useful  in  describing  the  early  phases  of  an  intense  explosion.  This  solu¬ 
tion,  is  of  a  similarity  nature,  i.e.  one  in  which  space  and  tiaie  are  simply 
related  and  in  a  manner  which  allows  the  partial  differential  equations  of 
motion  to  be  reduced  to  ordinary  differential  equations.  It  is  valid  for 
the  strong  shock  regime  of  poin.t  source  explosions  in  an  ideal  gas.  The 
most  widely  known  version  is  that  of  G.  I.  Taylor, however,  both  Sedov^^^ 

(3) 

in  Russia  and,  von  Neumann'  in  the  Ikiited  States,  at  about  the  same  time 
and  independently,  arrived  at  couparable  solutions.  These  were  all 
accomplished  during  World  War  II.  Aside  from  some  quite  different  ■ 
treatments  of  other  aspects  of  explosions  this  strong  shock,  point  soxirce, 
ideal  gas  solution  stood  without  significant  in^rovement  for  many  years 
as  the  best  description  of  the  early  phases  of  intense  explosions.  The 
Kirkwood -Brinkley  similarity  solution, invented  to  describe  b^eb 


P-1935 

2-16-60 

2 


explosive  blasts  also  origimted  during  World  War  II  as  did  a  •weak  shock, 

( 5) 

no  similarity  semiacoustical  treatment  by  Bethe.'  None  of  these  solutions, 
however,  was  of  sufficient  generality  to  meet  the  growing  requirements  for 
coB^lete  specification  of  blast  parameters,,  since  in  each  case  the  form  of 
the  solution  was  forced  by  the  similarity  assunption  or  (as  in  the  case  of  , 
Bethe's  solution)  the  limiting  physical  assumptions. 

Simi33r;!ty  solutions,  which  serve  to  reduce  the  equations  of  motion  to 
ordinai*y  differential  equations,  were  not  possible  outside  of  the  strong 
shock  reginie  except  for  very  special  circumstances,  von  Neumann  was 
among  the  first  to  recognize  the  necessity  of  resorting  to  numerical 
methods  for  the  more  general  solutions  of  the  blast  wave  problems.  Being 
also  a  pioneer  in  the  development  of  high  speed  electronic  computers,  he 
foresaw  the  usefulness  of  such  machines  in  pexfoming  truly  prodigious 
ninaerical  calculations.  His  continued  interest  in  both  machines  and 
blast  wave  problems  was  instrumental  to  the  developsient  of  a  number  of 
powerful  techniques.  In  any  methods  of  solution,  special  treatment  of 
shock  wave  propagation  is  required.  One  early  technique  was  to  fit 
separately  the  shock  position  1?y  means  of  the  governing  conservation 
laws  across  the  shock  front,  and  then  calculating  separately  the  effects 

(5) 

of  the  adiabatic  expansion  behind  the  shock. ITo  the  same  end  von 
Neumann  and  Richtmyer  introduced  artificial  viscosity,  the  effect  of 
\diich  was  to  spread  shock  discontinuities  in  such  a  way  that  they  could 

(7) 

be  treated  in  a  general  scheme.'  Later  a  partially  successful  method 

was  developed  by  lax  which  accooQ)lished  the  same  spreading  without  the 

(8) 

introduction  of  special  viscous  terms.'  '  Numerous  early  attenpts  at 
numerical  e<»putatlon0  proved  successful.  Notable  among  these  were 
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calculations  by  von  Neumann  and  Goldstein/^ ^  by  Wecken^^^  and  by  the 

( 7) 

present  author.  One  strong  impetus  for  the  developnent  of  methods  of 
conq?uting  general  hydrodynamic  shock  problems  was  in  the  developMnt  of 
nuclear  weapons,  and  much  of  the  pioneer  work  in  such  methods  originated 
in  the  los  Alamos  laboratory  of  the  ABC.  Considerable  interest  in  the 
effects  of  weapons,  both  nuclear  and  high  e^loslve,  pronqpted  further 
development  in  other  laboratories  and  institutions.  Attention  will  be 
focxised  in  the  following  section  on  one  such  method  of  approaching 
numerical  calculations  of  hydrodynamic  motions  in  conpressible  fluids.  ' 
The  intent  is  not  to  provide  the  audience  w;ith  an  ability  to  program 
and  run  shock  problems,  since  such  an  objective  should  require  a  more 
ambitious  course  and  would  lead  through  details  of  the  flow  diagraming  of 
the  program  logic  and  bookkeeping  methods.  A  portion  of  that  sort  of 
detail  may  be  found  in  Befs.  10  and  11.  Rather,  the  object  here  is  to 
acquaint  those  Interested  in  the  nature  of  numerical  methods  with  one 
such  method  in  sufficient  detail  to  allow  them  a  fair  evaluation  of  the 
effort  and/or  ease  of  formulating  such  techniques  and  of  using  such  a 
schfiane  to  solve  a  number  of  Interesting  problems. 
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11,  A  METHOD  07  HIMERICAL  CALCUUTION 

The  problem  presented  by  ejqplosions  lies  in  the  solution  of  a  set  of 
nonlinear,  partial  differential  equations  representing  the  conservation  of 
mass,  momentum  and  energy.  These  conservation  laws  may  be  stressed  matbe- 
matically  in  several  v&yB,  but  are  generally  formulated  in  terms  of  either 
Eulerian  or  Lagrangean  coordinates.  The  Eulerian  form  is  ein  eacpression  of 
the  conservation  laws  as  viewsd  from  coordinate  systems  fixed  in  space  and 
the  Lagrangean  form  is  an  ej^jression  of  the  same  conservations  in  terms  of 
a  fixed  set  of  masses  or  gas  particles,  so  that  a  solution  in  the  Eulerian 
case  represents  the  history  of  the  blast  wave  at  a  fixed  point  while  in  the 
Lagrangean  system  a  solution  describes  the  e:q>erience  of  a  particle  (or  a 
volume  of  gas)  as  it  moves  about.  Tbe  set  of  equations  below  is  one 
representation  of  these  conservation  laws  in  a  Lagrangean  system  for  a 
spherically  symmetric  explosion.  In  these  equations  the  lagrangean 
variable  or  mass  coordinate  is  characterized  by  the  mass  per  steradian, 

X 

i.e.,  m  «  PiRf/3  where  is  the  initial  density  of  the  gas  and  R  is  its 

X  Q  X  O 


initial  position. 

1  ^  1  hfP 

p*^  =  5  sr  ^ 

conservation  of  mass; 

(1) 

^  S  ^  » 

conseirvation  of  momentum  j 

(2) 

^  -  (P  +  Q)  ^  +  D  , 

conservation  of  energy^ 

(3) 

E  »  E(T,V),  P  =»  P(T,V)  1 

or 

P  »  P(E,V),  T  =  T(E,V)  I 

equations  of  state; 
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u  =  ^  ,  particle  velocity,  (5) 

In  these  equations  the  hynhols  are  defined  as  foUovs: 

V  -  specific  volume 
R  -  radial  distance 
u  -  particle  velocity 
P  -  press\ire 

Q  -  artificial  viscosity  pressure 
E  -  specific  internal  energy 
T  -  teB5>erature 

D  -  specific  energy  loss  (or  source)  rate 
m  -  aass/steradlan 
t  -  time 

In  solving  these  equations  for  a  specific  set  of  initial  and  boxmdary 
conditions  by  numerical  means,  the  usual  approach  is  to  arrange  a  set  of 
difference  equations  irtiich  approximate  these  differential  equations  over 
a  set  of  finite  mass  elements.  In  such  a  system,  then,  derivatives  become 
differences,  and  a  solution  is  accontplished  by  the  repetition  of  a  pre¬ 
scribed  series  of  numerical  steps  vhich  solve  the  difference  equations  for 
new  values  of  the  parameters  which  wo\ild  occur  after  a  small  time  increment. 

IDo  foimalize  the  differencing  procedure,  consider  the  transformation 
of  the  above  differential  equations  to  smother  set  in  'vdiich  the  independent 
variables  are  a  dimensionless  lagrangean  coordinate  j  and  a  dimensionless 
time  n,  and  in  which  the  masses  are  represented  by  Integer  values  of  J  and 
the  time  by  integer  values  of  n. 

The  Jacobian  for  the  lagrangean  coordinate  system  then  follows  from 
the  differential  equation  representing  conservation  of  mass. 
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1 

3V(j,0) 


a 


15 

R(j,o)  j 


1  a 

5V(j,n) 


R{j,n) 


3 


(6) 


The  other  conservation  equations  become 


^Sn^  ^  ^  (moaentum),  (7) 


^  -  (P  +  Q)  ^  +  D  (^)  (energy) , 


(8) 


and 


u  = 


aR  /at^ 


-1 


(9) 


In  these  equations  the  dex>eivdence  of  the  variables  on  J  and  n  has  been 
suppressed  for  sin^jlicity. 

Ihe  artificial  viscosity  may  be  variously  defined,  but  the  following 
is  a  usual  form: 


Q 


if 


Q  =  0 


if 


(10) 


in  which  the  value  of  C  determines  the  amount  of  spreadixig  of  shocks. 

When  C  is  unity,  the  spread  is  about  three  zones,  end’  when  C  is  six  the 
spread  is  around  six  zones,  the  constant  being  proportional  to  the  8q;uare 
of  the  spread. 

In  sone  unusual  clrcumsteinces  this  form  falls  to  provide  the  necessary 
shock  spreading.  Such  a  case  can  arise  in  some  problems  involving  Ijqploding 


shock  waves 
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A  more  rlgoroiis  foim  for  the  artificial  viscosity,  ii^ch  is  generally 
equivalent  to  the  form  above  (Eq.  (lO)),  is  the  following: 


« =  (|>'  ±  (i)"  (s)’" ' 


O  V 

Q  »  0  for  ^  ^  0 . 


The  formation  of  a  set  of  centered  difference  equations  is  now  a 
matter  of  assigning  integer  values  to  the  independent  variables  j  and  n 
and  either  defining  the  dei>endent  variables  at  these  points  or  at  midpoints 
(at  half  integer  points)  in  j  or  n.  It  is  conventional  to  esqpress  the  j 
and  n  values  associated  with  a  given  parameter  by  subscripts  and  super¬ 
scripts,  thus  the  radius  at  a  mass  point  J  and  at  a  time  t  after  n  time 
increiDents  At(n)  =  At^  would  appear  as  R^. 

The  equations  for  the  n-th  time  step  At”^  over  a  set  of  mass  elements 
follow. 


^  ^  ^  '  (la, 


for  which 


‘=“j  “  I  ^ 


At”  =  I  *  i  At“-i  , 


At”^  , 


(13) 
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j^+1  ^  j  2  2  ( 


1  n+4  n+4  ^ 

,  L  J  J-lJ ,_ 

^  v“  X 
0*5  J"2 


_  !!+■?  . 
for  Uj  ^  <  Uj^f  , 


<?t|  -  °  >  “"t!  - 


b“«  .  e“  (i  *  i  p“  .  *  -  v“  .)  f  D”i  .  (16) 

0“5  J"t  2  O-J  2  j-2  i'ls  J-? 

In  these  equations  it  is  warth  noting  that  the  quantities  are  either 
defined  or  averaged  in  such  a  uay  that  first  order  errors  in  At  or  Z!bi 
are  avoided,  i.e.,  the  difference  equations  are  centered.  The  actual 
progression  of  the  solution  can  better  he  illustrated  hy  reference  to  a 
space -time  grid,  as  helov,  showing  the  points  at  which  the  various  para- 


meters  are 

defined. 

n+1 

- R- 

t 

1 

n+j 

( 

- -  u  - 

1 

1 

n 

1 

- R- 

1 

n-j 

\ 

\ 

- u  - 

» 

1 

( 

j-i 

V,P,E - R - V,P,E 


- Q,D - u - Q,D-  - 


—  y,p,E- 


[ - V,P,E - 


__  - u - Q,D-- 


J+1 


For  illustration,  assmae  that  the  n-th  time  step  has  been  coinpleted,  and 
that  therefore  all  the  values  at  tijnes  less  than  n  are  known,  and  that  the 
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next  step  is  to  deteimine  the  valiies  at  ri+^  and  n+1  times.  The  difference 
equations  have  already  been  written  down  in  the  order  in  which  they  are 
computed.  The  first  equation  defines  a  new  velocity  at  time  n+5,  at  point 
in  terms  of  the  old  velocity  at  that  j -point  at  time  n-|-  and  the  radius 
at  that  j-point  at  time  n,  together  with  the  pressure  gradients  at  that 
point  and  at  time  n  as  represented  by  the  differences  in  the  pressures  at 
J  points  ^  a  unit  on  either  side  of  J.  From  this  new  velocity  a  new 
radius  is  coii5)uted  in  a  straightforward  manner,  using  the  definition  of 
the  particle  velocity.  These  new  radii,  in  turn,  are  used  to  coaqpute  a 
new  specific  voluime,  and  that  new  volume  together  with  the  new  values  of 
the  velocity  C6in  be  xised  to  cojxpute  a  new  airbificial  viscosity.  The 
artificial  viscosity  is  nearly  of  second  order  in  the  tjjne,  and  its  lack 
of  symmetry  in  the  velocity  equation  does  not  disturb  the  degree  of 
accuracy  of  these  otherwise  centered  difference  equations.  The  new 
volume,  together  with  the  new  artificial  viscosity,  can  then  be  used  in 
the  energy  equation  to  determine  a  new  energy  and  a  new  pressure.  If  the 
equation  of  state  is  not  siiaple,  it  is  then  necessary  to  arrive  at  the  new 
energy  and  pressure  by  some  iterative  method.  The  equation  of  state 
represents  a  second  relation  between  these  remaining  two  unknowns,  i.e., 
the  new  energy  and  the  new  pressure.  With  a  satisfactory  solution  for  all 
of  these  parameters  at  a  point  j  one  can  then  advance  to  the  next  J-point 
and  repeat  the  same  procedure. 

Of  course,  if  the  problem  involves  other  mechanisms  of  energy  transfer, 
or  sources  or  sinks  or  energy,  the  term  0  in  the  energy  equation  becomes 
useful.  If  this  term  is  to  represent  the  diffusion  of  energy,  either  by 
radiation  or  by  conduction,  then  ine  nature  of  the  problem  is  changed,  and 
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in  most  cases  is  made  much  more  conplicated.  If  for  instance,  ‘radiation 
diffusion  is  to  be  considered,  then  the  term  D  can  be  ejqjressed  as  a 
gradient  of  a  radiation  energy  flux* 


where 


“  ■  -  - 


(17) 


in  which  k{T,V)  represents  a  mean  opacity,  a  is  the  radiation  density 
constant  and  c  is  the  speed  of  light*  The  inclusion  of  such  terms  however, 
couples  the  difference  equations  for  energy  conservation  at  each  mass 
point,  so  that  one  then  is  faced  with  the  solution  of  something  like  two-j 
simultaneous  nonlinear  equations,  axid  the  iterative  procedure  must  become 
correspondingly  more  elegant. 

The  D  term  can  also  be  used  when  a  source  of  energy  is  required  as  an 
energy  rate  term  or,  when  energy  losses  occxir,  as  an  energy  sink*  For 
instance,  the  D  term  can  be  used  satisfactorily  in  detonation  wavs  calcula¬ 
tions*  A  nxsaber  of  forms  to  represent  the  bximing  of  high  explosives  or 
gases  in  a  detonation  front  are  possible;  one  fairly  crude  but  workable 
form  is  the  following: 


At 


n+^. 


CJ 


(18) 


s  m 


±k 


cj 


r 


in  which  is  the  Chapman-Jouget  energy,  i*e.,  the  energy  generated  in 
the  detonation,  "cj  is  the  detoi»tion  velocity,  and  S  is  the  spread  oi’ 
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thickness  of  the  detonation  front  in  keeping  with  the  shock  spreading 
characteristic  of  the  artificial  viscosity  technique. 

A  great  nany  cautions  resulting  fran  e^qperlence  and  from  the  details 
of  the  handling  of  the  numerical  procedure  might  well  be  included  in  a 
complete  description  of  a  numerical  method,  since  many  special  considera¬ 
tions  are  necessary  before  a  numerical  program  becMies  truly  workable. 

Some  of  these  factors  are  concerned  with  the  stability  of  the  numerical 
procedure,  attention  to  bovmdary  conditions  and  initial  conditions, 
techniques  for  expanding  the  grid  system  as  shocks  progress,  and  all  the 
many  probleias  of  input  and  output,  the  storing  of  information,  the 
stopping  and  restarting  of  the  problem,  and  internal  checks  on  the  accuracy 
of  the  computation. 

The  stability  criteria  which  must  be  considered  are  essential  limita¬ 
tions  on  the  size  of  the  increments  that  can  be  permitted.  The  Cpurant 
Condition  is  a  restriction  on  the  time  increment  (At)  such  that  local 
sounds  or  signals  from  one  mass  point  will  not  have  had  time  to  reach 
beyond  adjacent  mass  points  in  that  leisgth  of  time. 


n+1 


n 


\ 

\ 

> 

/ 

\ 

/ 

\ 

/ 

\ 

/ 

\ 

/ _ 

j-1  j  AR  j+i 


At 


At 


or 


At  <  —- 
c 
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For  an  ideal  gas  c  ■  7FV  and  generally  R  so  that  a  sufficient 

condition  for  stability  is 


At 


n+j 


VAbl  I 


(r^V7p/v/^„ 


(19) 


%diere  7  is  at  least  as  large  as  the  largest  possible  value  of  7  as 
determined  by  the  equation  of  state.  In  shock  fronts  (or  equivalently 
in  compression  regions)  the  artificial  viscosity  changes  the  nature  of 
the  differential  equations  from  that  of  a  wave  equation  to  that  of  a 
diffusion-type  equation,  and  another  stability  condition  must  be  included 
on  this  account.  An  approximate  derivation  of  this  stability  condition 
comes  from  the  momentum  equation  in  v^ich 


-  R‘ 


a 


T 


f 


and  from  the  conservation  of  mass. 


and 


R 


SR 


or 


i  £  1 

AR  V 


^  C*  (Au)^ 


or 


^  ai<o) 


for  diffusion  stability  then, 


At  <  ( 


AR 


C'/Au 


■)  (An)‘ 


Various  approximate 
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forms  of  this  inequality  can  be  used,  but  perhaps  the  siioplest  results  in 


2C 


^4  - 


<  1 


(20) 


\d3ere  C  is  the  constant  in  the  artificial  viscosity  (Eq.  (lO)). 

The  bovmdary  conditions  at  the  origin  in  cylindrical  or  spherical 
geometries  are  particvilarly  easily  satisfied  with  the  niaserical  scheme 
dealt  with  here,  since  it  is  only  necessary  to  note  that  the  velocity 
and  radi\i8  is  initially  zero  at  the  origin  emd  then  to  avoid  changing 
them  in  the  con^putation  loops  by  always  starting  a  cycle  with  the  concpute 
of  the  velocity  and  radixis  at  the  first  point  away  from  the  origin.  The 
artificial  viscosity  removes  the  necessity  of  any  complications  with 
floating  boundary  conditions  such  as  a  shock  front  might  otherwise  present, 
auad  the  other  boimdary  conditions  in  a  lagrangean  system  are  most  often 
straightforwardly  expressed  from  tlje  physical  restraints  of  the  problem. 

Fbr  instance,  if  the  calculation  is  of  an  explosion  in  air,  the  boundary 
condition  at  the  front  of  the  ejqplosion  may  easily  be  the  conditions  of 
the  ambient  air  in  undisturbed  zones  ahead  of  any  shock.  When  the  last 
zone  is  about  to  be  Influenced  by  the  shock,  e.g.,  its  temperature , or 
velocity  are  about  to  rise,  then  another  undisturbed  zone  can  be  added 
if  it  is  desirable  to  keep  a  limited  number  of  mass  zones  in  the  problem, 
then  a  scheme  can  be  arranged  whereby  when  a  zone  is  added  some  other 
pair  of  zones  is  combined  into  one  zone,  keeping  the  number  of  zones 
fixed.  Care  should  be  taken  in  such  an  exercise  to  observe  the  usual 
conservations  of  mass,  momentum  and  energy  by  altering  the  velocities 
and  internal  energies  of  the  resulting  coablned  zone  accordingly* 
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Stating  botmdary  conditions  for  moving  pistons  or  for  free  siirfaces 
bloving  off  into  a  vacuxoa  involves  no  more  serious  cobqplication  to  the 
code  than  to  alter  the  details  of  how  each  cycle  is  begun  or  ended*  For 
exajii5>le  in  the  blow-off  situation  each  cycle  is  stopped  at  the  last  zone, 
regardless  of  its  velocity  or  teniperature.  For  piston  motion,  on  the  other 
hand,  a  separate  prescription  must  be  provided  for  the  piston  velocity  or 
its  position,  but  then  the  usual  equations  are  followed,  for  each  mass 
point  =  f(t  or  n),  u^  =  g(t  or  n)»  !Ehe  specification  of  initial 
conditions  (or  botuadary  conditions,  for  that  matter)  is  a  mtter  of 
stating  the  physics  of  a  problem  in  terms  of  initial  values  of  a  set  of 
representative  variables  such  as  the  velocities,  temperatures  and  densities 
of  all  the  mass  eleiaents  at  some  Instant  of  time. 

Shells  or  layers  of  different  materials,  having  different  densities 
and  thermodynamic  properties  can  be  most  easily  managed  on  a  lagrangean 
system,  since  keeping  the  materials  separate  and  treating  each  differently 
only  involves  using  different  equations  of  state  according  to  which  j- 
values  (mass  elements)  the  coqputation  is  working  on. 

Numerical  techniques  of  the  form  described  here  possess  such 
generality  that  results  may  include  quite  unexpected  features  such  as 
unpredioted  shocks.  So  much  freedom  or  lack  of  forcing  in  the  form  or 
natirre  of  solutions  makes  the  method  at  once  very  powrful  and  yet  injects 
more  of  an  experimental  physics  flavor  into  the  procedure  in  contrast  to 
the  physical  clarity  of  analytic  results.  For  this  very  reason,  the 
processing  and  ultimate  analysis  of  numerical  results  is  a  vital  part  of 
this  approach  to  problem  solving  in  physics. 
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